import numpy as np
a = np.random.randn(10,10)
b= np.random.randn(10)
a @ b
array([-0.63269674, -1.32769635, -0.09029985, -3.52182299, -6.03363023,
        3.01976296, -2.75007156,  0.22389647, -1.55148195, -1.62650482])
b.T @ a.T
array([-0.63269674, -1.32769635, -0.09029985, -3.52182299, -6.03363023,
        3.01976296, -2.75007156,  0.22389647, -1.55148195, -1.62650482])
np.eye(10)
array([[1., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 1., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 1., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 1., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 1., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 1., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 1., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 1., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 1., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 1.]])
np.ones((10,10))
array([[1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1., 1., 1., 1.]])
np.linalg.inv(a) @ a
array([[ 1.00000000e+00, -3.65595152e-16, -1.37241160e-16,
        -1.67938840e-15,  9.33527437e-16,  1.44199349e-15,
         1.15089786e-16, -1.55696214e-15,  1.11022302e-15,
        -2.63677968e-16],
       [ 1.26550255e-15,  1.00000000e+00, -2.68494358e-16,
         1.42303383e-15, -9.79693103e-16, -4.51521472e-15,
         2.63044668e-15,  1.74407151e-15, -2.66453526e-15,
         4.44089210e-16],
       [-1.64423122e-15, -1.66173972e-16,  1.00000000e+00,
         2.75288470e-15,  1.74102465e-16,  1.54544175e-15,
        -2.79739957e-15, -2.04740762e-15,  0.00000000e+00,
         8.88178420e-16],
       [-3.51858931e-16,  6.23676926e-17,  2.66458886e-16,
         1.00000000e+00,  9.38263196e-17,  2.61301662e-16,
         2.45064178e-16, -1.53357911e-16,  4.44089210e-16,
        -4.44089210e-16],
       [ 7.84703105e-16, -1.04920005e-15,  4.13020606e-16,
        -2.12483741e-15,  1.00000000e+00, -4.53483867e-15,
         4.18017106e-15,  1.71820802e-15,  0.00000000e+00,
         0.00000000e+00],
       [ 9.37535921e-16,  3.43043467e-16,  9.65347805e-17,
         6.68706576e-16, -4.20016708e-16,  1.00000000e+00,
        -1.09267442e-15,  1.01169627e-16,  0.00000000e+00,
         4.44089210e-16],
       [-8.71025446e-16,  1.09123436e-15, -7.64517688e-16,
         9.73897302e-16, -6.40358430e-16,  3.68586912e-15,
         1.00000000e+00, -1.54120216e-15,  0.00000000e+00,
         8.88178420e-16],
       [ 2.86959469e-15,  2.45431198e-16,  5.06260913e-16,
        -5.81240373e-16, -5.42632247e-16, -2.38411883e-15,
         1.09159912e-15,  1.00000000e+00, -8.88178420e-16,
        -4.44089210e-16],
       [-3.57782462e-17, -5.09336718e-16,  7.90198623e-16,
        -2.53964484e-15, -7.71687051e-17, -1.71387168e-15,
         1.28935306e-15,  1.41638585e-15,  1.00000000e+00,
        -8.88178420e-16],
       [-5.10166112e-16,  3.36418889e-16,  3.54329540e-16,
        -1.34491999e-15, -3.29604245e-16, -7.44603616e-16,
         1.18536447e-15,  5.30313254e-16,  4.44089210e-16,
         1.00000000e+00]])
np.linalg.svd(a).Vh
array([[ 0.28383114,  0.13031739,  0.08488282,  0.26915419, -0.05833712,
        -0.55745972,  0.04622785,  0.52583554, -0.47933131, -0.01565541],
       [-0.7299212 ,  0.04346431, -0.15481374,  0.1199476 ,  0.28639583,
        -0.35620311,  0.24101023, -0.22640736, -0.21810744, -0.24728702],
       [ 0.07621821,  0.43546203,  0.2088927 ,  0.28230896, -0.17741055,
        -0.17934521,  0.42879952, -0.01233035,  0.62390845, -0.21049244],
       [-0.39305518,  0.11448225,  0.09008748,  0.41564326, -0.41297868,
         0.10835773, -0.0692628 , -0.00268073, -0.05982432,  0.67885942],
       [ 0.06933558,  0.34860056, -0.39022658,  0.26430987, -0.41367073,
         0.19057494, -0.35422431, -0.22534478, -0.22225331, -0.46737544],
       [-0.05118823, -0.68217798, -0.03157757,  0.58681909,  0.00875722,
         0.07532819, -0.10790011,  0.20278167,  0.24018863, -0.26543311],
       [ 0.31944117, -0.2503526 , -0.62925149,  0.04300488, -0.15611409,
        -0.20160406,  0.45924193, -0.30071194, -0.00703338,  0.26660451],
       [ 0.25406731,  0.17492689,  0.15182522,  0.44334045,  0.50247181,
         0.45816405,  0.28794845, -0.19810239, -0.31166484,  0.05968817],
       [ 0.22725739, -0.10877031,  0.36247981,  0.15032739,  0.06499115,
        -0.45411298, -0.38287566, -0.64660084,  0.0248702 ,  0.08237576],
       [ 0.02620891,  0.29393688, -0.46195436,  0.15967492,  0.51017133,
        -0.14029136, -0.41663193,  0.18003309,  0.35346832,  0.25124871]])
np.diagonal?
Signature:       np.diagonal(a, offset=0, axis1=0, axis2=1)
Call signature:  np.diagonal(*args, **kwargs)
Type:            _ArrayFunctionDispatcher
String form:     <function diagonal at 0x7b3729c832e0>
File:            ~/miniforge3/envs/stan/lib/python3.11/site-packages/numpy/core/fromnumeric.py
Docstring:      
Return specified diagonals.
If `a` is 2-D, returns the diagonal of `a` with the given offset,
i.e., the collection of elements of the form ``a[i, i+offset]``.  If
`a` has more than two dimensions, then the axes specified by `axis1`
and `axis2` are used to determine the 2-D sub-array whose diagonal is
returned.  The shape of the resulting array can be determined by
removing `axis1` and `axis2` and appending an index to the right equal
to the size of the resulting diagonals.
In versions of NumPy prior to 1.7, this function always returned a new,
independent array containing a copy of the values in the diagonal.
In NumPy 1.7 and 1.8, it continues to return a copy of the diagonal,
but depending on this fact is deprecated. Writing to the resulting
array continues to work as it used to, but a FutureWarning is issued.
Starting in NumPy 1.9 it returns a read-only view on the original array.
Attempting to write to the resulting array will produce an error.
In some future release, it will return a read/write view and writing to
the returned array will alter your original array.  The returned array
will have the same type as the input array.
If you don't write to the array returned by this function, then you can
just ignore all of the above.
If you depend on the current behavior, then we suggest copying the
returned array explicitly, i.e., use ``np.diagonal(a).copy()`` instead
of just ``np.diagonal(a)``. This will work with both past and future
versions of NumPy.
Parameters
----------
a : array_like
    Array from which the diagonals are taken.
offset : int, optional
    Offset of the diagonal from the main diagonal.  Can be positive or
    negative.  Defaults to main diagonal (0).
axis1 : int, optional
    Axis to be used as the first axis of the 2-D sub-arrays from which
    the diagonals should be taken.  Defaults to first axis (0).
axis2 : int, optional
    Axis to be used as the second axis of the 2-D sub-arrays from
    which the diagonals should be taken. Defaults to second axis (1).
Returns
-------
array_of_diagonals : ndarray
    If `a` is 2-D, then a 1-D array containing the diagonal and of the
    same type as `a` is returned unless `a` is a `matrix`, in which case
    a 1-D array rather than a (2-D) `matrix` is returned in order to
    maintain backward compatibility.
    If ``a.ndim > 2``, then the dimensions specified by `axis1` and `axis2`
    are removed, and a new axis inserted at the end corresponding to the
    diagonal.
Raises
------
ValueError
    If the dimension of `a` is less than 2.
See Also
--------
diag : MATLAB work-a-like for 1-D and 2-D arrays.
diagflat : Create diagonal arrays.
trace : Sum along diagonals.
Examples
--------
>>> a = np.arange(4).reshape(2,2)
>>> a
array([[0, 1],
       [2, 3]])
>>> a.diagonal()
array([0, 3])
>>> a.diagonal(1)
array([1])
A 3-D example:
>>> a = np.arange(8).reshape(2,2,2); a
array([[[0, 1],
        [2, 3]],
       [[4, 5],
        [6, 7]]])
>>> a.diagonal(0,  # Main diagonals of two arrays created by skipping
...            0,  # across the outer(left)-most axis last and
...            1)  # the "middle" (row) axis first.
array([[0, 6],
       [1, 7]])
The sub-arrays whose main diagonals we just obtained; note that each
corresponds to fixing the right-most (column) axis, and that the
diagonals are "packed" in rows.
>>> a[:,:,0]  # main diagonal is [0 6]
array([[0, 2],
       [4, 6]])
>>> a[:,:,1]  # main diagonal is [1 7]
array([[1, 3],
       [5, 7]])
The anti-diagonal can be obtained by reversing the order of elements
using either `numpy.flipud` or `numpy.fliplr`.
>>> a = np.arange(9).reshape(3, 3)
>>> a
array([[0, 1, 2],
       [3, 4, 5],
       [6, 7, 8]])
>>> np.fliplr(a).diagonal()  # Horizontal flip
array([2, 4, 6])
>>> np.flipud(a).diagonal()  # Vertical flip
array([6, 4, 2])
Note that the order in which the diagonal is retrieved varies depending
on the flip function.
Class docstring:
Class to wrap functions with checks for __array_function__ overrides.
All arguments are required, and can only be passed by position.
Parameters
----------
dispatcher : function or None
    The dispatcher function that returns a single sequence-like object
    of all arguments relevant.  It must have the same signature (except
    the default values) as the actual implementation.
    If ``None``, this is a ``like=`` dispatcher and the
    ``_ArrayFunctionDispatcher`` must be called with ``like`` as the
    first (additional and positional) argument.
implementation : function
    Function that implements the operation on NumPy arrays without
    overrides.  Arguments passed calling the ``_ArrayFunctionDispatcher``
    will be forwarded to this (and the ``dispatcher``) as if using
    ``*args, **kwargs``.
Attributes
----------
_implementation : function
    The original implementation passed in.
np.diag([1,2,3])
array([[1, 0, 0],
       [0, 2, 0],
       [0, 0, 3]])
b = np.random.randn(100,2)
np.linalg.svd(b)
SVDResult(U=array([[ 3.18048731e-05,  6.74648048e-02, -6.41579788e-02, ...,
        -4.70143275e-02, -1.17012062e-02,  5.72707204e-02],
       [-1.67854420e-01,  1.70642880e-01, -1.21850722e-01, ...,
        -1.51615802e-01,  2.00743368e-02, -1.09802273e-01],
       [-7.00863302e-02,  1.18377452e-01,  9.84604342e-01, ...,
        -1.72681740e-02,  1.25486957e-03, -7.25125082e-03],
       ...,
       [-1.04494452e-01,  1.19440600e-01, -1.73568100e-02, ...,
         9.79733632e-01,  1.95664282e-03, -1.09753255e-02],
       [ 2.32789965e-02, -5.98960947e-04,  1.31502367e-03, ...,
         2.01688869e-03,  9.99525031e-01,  2.51988849e-03],
       [-1.23784387e-01,  8.99047891e-03, -7.56209658e-03, ...,
        -1.12822312e-02,  2.51689274e-03,  9.86617378e-01]]), S=array([9.7586975 , 8.78873925]), Vh=array([[-0.46515128, -0.8852312 ],
       [-0.8852312 ,  0.46515128]]))

\[ A = U\Sigma V^T\] \[ (A^TA)^{-1} A = (V\Sigma^T U^T U \Sigma^T V^T)^{-1} U\Sigma V^T \] \[A^+ = (V \Sigma^T\Sigma V^T)^{-1} V\Sigma^T U^T \] \[ A^+ = V (\Sigma^T\Sigma)^{-1}\Sigma^T U^T \] \[ A^+ = V\Sigma^+ U^T \]

U = np.linalg.svd(b)[0]
U.T @ U
array([[ 1.00000000e+00, -1.24938610e-16,  5.63615904e-17, ...,
         2.77555756e-17, -5.20417043e-18,  0.00000000e+00],
       [-1.24938610e-16,  1.00000000e+00, -3.15292394e-17, ...,
        -3.46944695e-17,  6.50521303e-18, -9.54097912e-18],
       [ 5.63615904e-17, -3.15292394e-17,  1.00000000e+00, ...,
        -3.46944695e-18, -1.30104261e-18,  5.20417043e-18],
       ...,
       [ 2.77555756e-17, -3.46944695e-17, -3.46944695e-18, ...,
         1.00000000e+00,  4.82754810e-19, -4.24892617e-18],
       [-5.20417043e-18,  6.50521303e-18, -1.30104261e-18, ...,
         4.82754810e-19,  1.00000000e+00,  3.18478586e-19],
       [ 0.00000000e+00, -9.54097912e-18,  5.20417043e-18, ...,
        -4.24892617e-18,  3.18478586e-19,  1.00000000e+00]])
$$ V(\Sigma^T\Sigma)^{-1}\Sigma^TU^T b = x $$
$$ (\Sigma^T\Sigma)^{-1}\Sigma^T U^T b  = V^T x$$
$$ (\Sigma^T U^T b) = (\Sigma^T\Sigma) (V^T x)$$

`V @ np.linalg.solve(S^T @ S, U.T @ b)